# a script that extracts data from the structure arrays, compiled by the "KORUS_rc_chart" script, and creates histogram plots for selected industries or GDP (only for the results of additive decomposition, and only for NET results)
# NOTE: "KORUS_rc_chart" script must be run first, so that all structures and years y0, y1 are available for a particular yearly interval
# version: Dec 2020
# ENSURE that Octave can find the script, use the "addpath" command: addpath ("/[insert the path to your Octave folder]/Octave")

# ask for input identifiers

# ask whether the histogram should be created for Korea or the U.S.
i = menu ("choose the country for the histogram:", "KOR", "USA");
if (i == 1)
    rs = ["r36_s19"];
    sr = ["r19_s36"];
    sname = ["KOR"];
elseif (i == 2)
    rs = ["r19_s36"];
    sr = ["r36_s19"];
    sname = ["USA"];
endif

# identify industry or GDP for which the histogram is created
k = menu ("choose the industry or GDP for the histogram:", "D01T03 Agriculture, forestry and fishing", "D05T06 Mining and extraction of energy producing products", "D07T08 Mining and quarrying of non-energy producing products", "D09 Mining support service activities", "D10T12 Food products, beverages and tobacco", "D13T15 Textiles, wearing apparel, leather and related products", "D16 Wood and products of wood and cork", "D17T18 Paper products and printing", "D19 Coke and refined petroleum products", "D20T21 Chemicals and pharmaceutical products", "D22 Rubber and plastic products", "D23 Other non-metallic mineral products", "D24 Basic metals", "D25 Fabricated metal products", "D26 Computer, electronic and optical products", "D27 Electrical equipment", "D28 Machinery and equipment, nec", "D29 Motor vehicles, trailers and semi-trailers", "D30 Other transport equipment", "D31T33 Other manufacturing; repair and installation of machinery and equipment", "D35T39 Electricity, gas, water supply, sewerage, waste and remediation services", "D41T43 Construction", "D45T47 Wholesale and retail trade; repair of motor vehicles", "D49T53 Transportation and storage", "D55T56 Accommodation and food services", "D58T60 Publishing, audiovisual and broadcasting activities", "D61 Telecommunications", "D62T63 IT and other information services", "D64T66 Financial and insurance activities", "D68 Real estate activities", "D69T82 Other business sector services", "D84 Public admin. and defence; compulsory social security", "D85 Education", "D86T88 Human health and social work", "D90T96 Arts, entertainment, recreation and other service activities", "GDP");
# NOTE: the list of industries above excludes "D97T98 Private households with employed persons" because it has usually 0 values and includes "GDP" at No.36

# ask whether to include total demand or only intermediate or final demand
input_dem = menu ("choose the composition of demand:", "total = intermediate + final", "intermediate", "final");
if (input_dem == 1)
    dem = ["AF"];
elseif (input_dem == 2)
    dem = ["A"];
elseif (input_dem == 3)
    dem = ["F"];
endif

# ask which effect should be decsribed by the histogram
input_eff = menu ("choose effect:", "trade creation (TC)", "trade diversion in (TDin)", "trade diversion out (TDout)", "trade contraction (TR)");
if (input_eff == 1)
    eff = ["TC"];
elseif (input_eff == 2)
    eff = ["TDin"];
elseif (input_eff == 3)
    eff = ["TDout"];
elseif (input_eff == 4)
    eff = ["TR"];
endif

# for GDP, k = 36
if (k == 36)
# transform the target structure array into a matrix (additive SDA results)
# for the target country as home
    z = 1;
    for [val, key] = ALLdy.(eff).(dem).(rs)
        Myrs(:, z++) = val/1000;
    endfor
# for the target country as partner
    z = 1;
    for [val, key] = ALLdy.(eff).(dem).(sr)
        Mysr(:, z++) = val/1000;
    endfor
# add up to obtain the net results
    my = Myrs(1,:) + Mysr(2,:);
# create histogram with 20 bins
    hist(my,20,"facecolor", [0.8, 0.8, 0.8], "edgecolor", "black")
    set(1,"papersize",[6 4.5])
    set(1,"paperposition",[0 0 6 4.5])
    set(gca(), "fontsize", 16)
    ylabel("count", "fontsize", 16)
    xlabel("change in GDP, $ billion", "fontsize", 16)
# calculate the average (the baseline estimate)
    m = mean(my);
# set the x coordinates for the beginning and end of the line - average
    x = [m, m];
# get the histogram counts and centers of bins
    [nn, xx] = hist(my,20);
# identify the bin that contains the average
    bin = lookup(xx - m, min(abs(xx - m)));
# set the y coordinates for the beginning and end of the line - average
    y = [0, nn(bin)];
# draw the line that denotes the baseline estimate, fitting the relevant histogram bin
    line(x,y,"color","blue", "linewidth",2);
# set the x and y coordinates for the beginning and end of the line - alternative estimate (from the simplfied, bi-proportional allocation of market share changes)
# for the target country as home
    bpyrs = BPdy.(eff).(dem).(rs);
# for the target country as partner
    bpysr = BPdy.(eff).(dem).(sr);
# add up to obtain the net results
    bpy = (bpyrs(1) + bpysr(2))/1000;
# set the x and y coordinates
    x = [bpy, bpy];
    bin = lookup(xx - bpy, min(abs(xx - bpy)));
    y = [0, nn(bin)];
# draw the line that denotes the alternative estimate, fitting the relevant histogram bin
    line(x,y,"color","red", "linewidth",2);
# save chart
    filename = ["/[insert the path to the folder with the charts, a subfolder for GDP is recommended]/charts (octave)_dy/Fig_" y0 "_" y1 "_dy_" eff "_" dem "_" sname ".pdf"];
    print(filename, "-dpdfcrop")
    clf
# clear temporary elements
    clear -x ALL* BP* y0 y1

# for a selected industry, k = [1,...,35]
elseif

# transform the target structure array into 3-d tensors (value added - additive)
# for the target country as home
    z = 1;
    for [val, key] = ALLdv.(eff).(dem).(rs)
        Tvrs(:, :, z++) = val;
    endfor
# for the target country as partner
    z = 1;
    for [val, key] = ALLdv.(eff).(dem).(sr)
        Tvsr(:, :, z++) = val;
    endfor
# rotate the tensors Tv[rs,sr] so that the first dimension is for country (i=1 for home country and i=2 for partner country), the second dimension is for the Monte Carlo runs (j=1,...,400) and the third dimension is for industry (k=1,...,36)
    Tvrs = flip(rotdim(Tvrs, -1, [2,3]),3);
    Tvsr = flip(rotdim(Tvsr, -1, [2,3]),3);
# extract the relevant data and add up to obtain the net results
    ty = Tvrs(1,:,k) + Tvsr(2,:,k);
# create a histogram with 20 bins
    hist(ty,20,"facecolor", [0.8, 0.8, 0.8], "edgecolor", "black")
    set(1,"papersize",[6 4.5])
    set(1,"paperposition",[0 0 6 4.5])
    set(gca(), "fontsize", 16)
    ylabel("count", "fontsize", 16)
    xlabel("change in value added, $ million", "fontsize", 16)
# calculate the average (the baseline estimate)
    m = mean(ty);
# set the x coordinates for the beginning and end of the line
    x = [m, m];
# get the histogram counts and centers of bins
    [nn, xx] = hist(ty,20);
# identify the bin that contains the average
    bin = lookup(xx - m, min(abs(xx - m)));
# set the y coordinates for the beginning and end of the line
    y = [0, nn(bin)];
# draw the line that denotes the baseline estimate, fitting the relevant histogram bin
    line(x,y,"color","blue", "linewidth",2);
# set the x and y coordinates for the beginning and end of the line - alternative estimate (from the simplfied, bi-proportional allocation of market share changes)
# for the target country as home
    BPvrs = BPdv.(eff).(dem).(rs);
# for the target country as partner
    BPvsr = BPdv.(eff).(dem).(sr);
# add up to obtain the net results
    bpv = BPvrs(1,k) + BPvsr(2,k);
# set the x and y coordinates
    x = [bpv, bpv];
    bin = lookup(xx - bpv, min(abs(xx - bpv)));
    y = [0, nn(bin)];
# draw the line that denotes the alternative estimate, fitting the relevant histogram bin
    line(x,y,"color","red", "linewidth",2);
# save chart
    filename = ["/[insert the path to the folder with the charts, a subfolder for each industry is recommended]/charts (octave)_dv/Fig_" y0 "_" y1 "_dv_ind" num2str(k) "_" eff "_" dem "_" sname ".pdf"];
    print(filename, "-dpdfcrop")
    clf

# repeat the above steps for employment - additive only, transform the target structure arrays into a 3-d tensor
# for the target country as home
    z = 1;
    for [val, key] = ALLde.(eff).(dem).(rs)
        Ters(:, :, z++) = val;
    endfor
# for the target country as home
    z = 1;
    for [val, key] = ALLde.(eff).(dem).(sr)
        Tesr(:, :, z++) = val;
    endfor
# rotate the tensors Te so that the first dimension is for country (i=1 for home country and i=2 for partner country), the second dimension is for the Monte Carlo runs (j=1,...,400) and the third dimension is for industry (k=1,...,36)
    Ters = flip(rotdim(Ters, -1, [2,3]),3);
    Tesr = flip(rotdim(Tesr, -1, [2,3]),3);
# extract the relevant data and add up to obtain the net results
    te = Ters(1,:,k) + Tesr(2,:,k);
# create a histogram with 20 bins
    hist(te,20,"facecolor", [0.8, 0.8, 0.8], "edgecolor", "black")
    set(gca(), "fontsize", 16)
    ylabel("count", "fontsize", 16)
    xlabel("change in employment, thousand persons/jobs", "fontsize", 16)
# calculate the average (the baseline estimate)
    m = mean(te);
# set the x coordinates for the beginning and end of the line
    x = [m, m];
# get the histogram counts and centers of bins
    [nn, xx] = hist(te,20);
# identify the bin that contains the average
    bin = lookup(xx - m, min(abs(xx - m)));
# set the y coordinates for the beginning and end of the line
    y = [0, nn(bin)];
# draw the line that denotes the baseline estimate, fitting the relevant histogram bin
    line(x,y,"color","blue", "linewidth",2);
# set the x and y coordinates for the beginning and end of the line - alternative estimate (from the simplfied, bi-proportional allocation of market share changes)
# for the target country as home
    BPers = BPde.(eff).(dem).(rs);
# for the target country as partner
    BPesr = BPde.(eff).(dem).(sr);
# add up to obtain the net results
    bpe = BPers(1,k) + BPesr(2,k);
# set the x and y coordinates
    x = [bpe, bpe];
    bin = lookup(xx - bpe, min(abs(xx - bpe)));
    y = [0, nn(bin)];
# draw the line that denotes the alternative estimate, fitting the relevant histogram bin
    line(x,y,"color","red", "linewidth",2);
# save chart
    filename = ["/[insert the path to the folder with the charts, a subfolder for each industry is recommended]/charts (octave)_de/Fig_" y0 "_" y1 "_de_ind" num2str(k) "_" eff "_" dem "_" sname ".pdf"];
    print(filename, "-dpdfcrop")
    clf
# clear temporary elements
    clear -x ALL* BP* y0 y1

endif

# the structures are still available, manually clear if not required any longer

# END of script
